Głównym problemem, który my musimy zrozumieć, a protony (ogólnie: naładowane dodatnio, odpychające się elektrostatycznie jądra atomowe) pokonać, jest bariera kulombowska. Otóż na pierwszy rzut oka wydaje się, że energie kinetyczne ruchów termicznych są o rzędy wielkości za małe, aby protony w ogóle mogły się do siebie zbliżyć, nie mówiąc nawet o zajściu reakcji jądrowych.  Typowy potencjał oddziaływania wygląda następująco.

Plot[{UnitStep[r-1/4]/r +-2 UnitStep[1/4-r],Labeled[2,"E"]},{r,0,5},PlotStyle->{Red,{Black,Dashed}},PlotRange->All,ExclusionsStyle->Red]

Klasycznie, cząstka o energii E musi odbić się od bariery kulombowskiej  i nigdy nie przedostanie się do obszaru silnych oddziaływań jądrowych reprezentowanych przez studnię potencjału dla . Kluczowe to jest zjawisko tunelowania kwantowego. Obliczenie współczynnika przejścia przez barierę to klasyczne zagadnienie z mechaniki kwantowej. Zrobimy to korzystając z przybliżenia WKB:

We wzorze powyżej,  to promień jądra, natomiast  to tzw. klasyczny punkt powrotu czyli miejsce, gdzie cząstka odbiłaby się gdyby nie było efektów kwantowych. Wykonujemy całkę.

Integrate[Sqrt[kk/r-En],{r,r0,kk/En},Assumptions->r0>0&&En<1/r0&&En>0&&kk>0&&kk/r0>En]

Interesuje nas przypadek, gdy bariera jest bardzo wysoka, tj. energia E mała. Rozwijamy w szereg potęgowy.

Series[%,{En,0,5}]//Normal//First

Wstawiamy wzór na prawo Coulomba w układzie SI.

P=Exp[-2/(h/(2Pi))Sqrt[2m]%]/.kk->qe^2 Za Zb/(4Pi (Subscript[ε, 0]))

P/.En->1/2m v^2 (* Proszę pamiętać, że m jest tu masą zredukowaną ! *)

PowerExpand[%]

Współczynnik tunelowania P zapisujemy często jako:

,

gdzie  nazywamy parametrem Sommerfelda. Wzór wygląda szczególnie ładnie zapisany za pomocą bezwymiarowej stałej struktury subtelnej .

η=PowerExpand[Log[%]/(-2Pi)/.qe^2->α(4Pi (Subscript[ε, 0]) h/(2Pi)c)]

Inny sposób, szczególnie wygodny w astrofizyce, to użycie energii Gamova, :

 

gdzie:

EG=(Sqrt[En]Log[P])^2//PowerExpand

Ładnie wygląda energia Gamova wyrażona w jednostkach :

EG/(m c^2) /.qe->Sqrt[α(4Pi (Subscript[ε, 0]) h/(2Pi)c)]

Sqrt[%]//PowerExpand